In [ ]:
# Bibliotecas
import os # para manipulação de arquivos
import numpy as np # para operações numéricas
import librosa # para processamento de áudio
import matplotlib.pyplot as plt # para criar gráficos e visualizações
import IPython.display as ipd # para reproduzir áudio diretamente no Jupyter Notebook
from scipy.signal import stft, istft # short term fourier transform
import pandas as pd # para criar tabelas
import seaborn as sns
1. Carregamento e Pré-processamento do Áudio¶
In [15]:
# Constantes otimizadas para análise de áudio
taxa_amostragem = 22050 # 22050 Hz é um valor de referência
nperseg = 1024 # valor de referência
# Carregamento de um audio de exemplo
audio, sr = librosa.load("audio/Confutatis.wav", sr=taxa_amostragem)
ipd.Audio(data=audio, rate=taxa_amostragem)
Out[15]:
2. Cálculo do Espectrograma (STFT)¶
In [16]:
# Calcular STFT
frequencias, tempos, matriz_stft = stft(audio, fs=taxa_amostragem, nperseg=nperseg)
espectrograma_db = 20 * np.log10(np.abs(np.abs(matriz_stft)) + 1e-10) # conversão para decibels
fase = np.angle(matriz_stft) # necessário para a reconstrução do áudio
3. Decomposição SVD e Análise Espectral¶
In [17]:
# Função para desenhar espectrograma
def plotar_espectrograma(S_db, sr, title="", ax=None):
img = librosa.display.specshow(S_db,
sr=sr,
hop_length=512,
x_axis='time',
y_axis='log',
cmap='inferno',
vmin=-80, vmax=0,
ax=ax)
cbar = plt.colorbar(img, ax=ax, format='%+2.0f dB')
cbar.set_label('Magnitude (dB)')
if ax:
ax.set_title("Espectograma"+" {}".format(title))
ax.set_ylabel('Frequência [Hz]')
ax.set_xlabel('Tempo [s]')
else:
plt.title("Espectograma"+" {}".format(title))
plt.ylabel('Frequência [Hz]')
plt.xlabel('Tempo [s]')
# Criar subplots lado a lado
plt.figure(figsize=(10, 4))
# Plotar espectrograma original
plotar_espectrograma(espectrograma_db, taxa_amostragem)
# Ajustar layout
plt.tight_layout()
plt.show()
In [18]:
def decomposicaoSVD(espectrograma, k=None):
U, S, Vt = np.linalg.svd(espectrograma, full_matrices=False)
U = U[:, :k]
S = np.diag(S[:k])
Vt = Vt[:k, :]
return U, S, Vt
In [64]:
# Decomposição SVD do espectrograma
_, matriz_val_singulares, _ = decomposicaoSVD(espectrograma=espectrograma_db)
# Extrair os valores singulares
valores_singulares = np.diag(matriz_val_singulares)
# Calcular estatísticas
media = np.mean(valores_singulares)
desvio = np.std(valores_singulares)
valor_max = np.max(valores_singulares)
valor_min = np.min(valores_singulares)
mediana = np.quantile(valores_singulares,0.5)
quartil_1 = np.quantile(valores_singulares,0.25)
In [69]:
# Decomposição SVD do espectrograma
_, matriz_val_singulares, _ = decomposicaoSVD(espectrograma=espectrograma_db)
# Extrair os valores singulares
valores_singulares = np.diag(matriz_val_singulares)
# Calcular estatísticas
media = np.mean(valores_singulares)
desvio = np.std(valores_singulares)
valor_max = np.max(valores_singulares)
valor_min = np.min(valores_singulares)
mediana = np.quantile(valores_singulares,0.5)
quartil_1 = np.quantile(valores_singulares,0.25)
In [220]:
# Cálculo da proporção acumulada
proporcao_acumulada = np.cumsum(valores_singulares) / np.sum(valores_singulares)
# Encontrar o menor k tal que a proporção acumulada ≥ 95%
k_95 = np.argmax(proporcao_acumulada >= 0.95)
# Plot
plt.figure(figsize=(10, 4))
plt.plot(proporcao_acumulada, color='navy', linewidth=2, label='Energia acumulada')
# Linha vertical no k_95
plt.xticks(list(plt.xticks()[0]) + [k_95])
# Linha horizontal em 95%
plt.axhline(0.95, color='red', linestyle='--', linewidth=1.2, label='95% da energia')
# Rótulos e ajustes
plt.xlabel('Número de valores singulares $k$')
plt.ylabel('Proporção acumulada normalizada')
plt.title('Proporção acumulada dos valores singulares')
plt.grid(True, linestyle='--', alpha=0.6)
plt.legend()
plt.tight_layout()
plt.show()
In [235]:
# Paleta com alto contraste e baixa saturação
cor_media = '#5e81ac' # Azul acinzentado
cor_mediana = '#a3be8c' # Verde oliva suave
cor_quartil = '#d08770' # Laranja queimado
cor_min_max = '#b48ead' # Roxo malva apagado
# Plot
plt.figure(figsize=(10, 4))
sns.histplot(valores_singulares, kde=True, log_scale=True, color='#d8dee9', edgecolor='gray')
# Linhas verticais com legenda
plt.axvline(media, color=cor_media, linestyle='--', linewidth=2, alpha=0.9, label=f'Média (μ = {media:.2f})')
plt.axvline(mediana, color=cor_mediana, linestyle='-', linewidth=2, alpha=0.9, label=f'Mediana = {mediana:.2f}')
plt.axvline(quartil_1, color=cor_quartil, linestyle=':', linewidth=2, alpha=0.9, label=f'1º Quartil = {quartil_1:.2f}')
plt.axvline(valor_max, color=cor_min_max, linestyle='-.', linewidth=1.5, alpha=0.4, label=f'Máx = {valor_max:.2f}')
plt.axvline(valor_min, color=cor_min_max, linestyle='-.', linewidth=1.5, alpha=0.4, label=f'Mín = {valor_min:.2f}')
# Configurações visuais
plt.xlabel("Valores Singulares")
plt.ylabel("Frequência (escala log)")
plt.title("Distribuição dos Valores Singulares com Estatísticas")
plt.grid(True, linestyle='--', alpha=0.4)
plt.legend(loc='upper right')
plt.tight_layout()
plt.show()
In [166]:
# Criar o gráfico
plt.figure(figsize=(12, 4))
plt.plot(valores_singulares,".", color='#457b9d', linewidth=2, label="Espectro de Valores Singulares")
plt.grid(True)
# Linhas verticais de truncamento
locais_truncamento = [10,25,50,100,150,200,300,400]
for linha in locais_truncamento:
plt.axvline(x=linha, color='red', linestyle='--', linewidth=0.75)
# Configurar os ticks do eixo x
xticks_custom = np.unique(np.concatenate((np.arange(0, len(valores_singulares), 100), locais_truncamento, [len(valores_singulares)])))
yticks_custom = np.unique(np.concatenate((np.arange(0, valor_max, 10000), [valor_max])))
plt.xticks(xticks_custom, rotation=45)
plt.yticks(yticks_custom)
# Legenda e títulos
plt.yscale('log')
plt.legend(['Valores Singulares', 'Locais de Truncamento'], loc='upper right')
plt.title(f"Valores Singulares (n={len(valores_singulares)})")
plt.xlabel('Quantidade de Valores Singulares')
plt.tight_layout()
plt.show()
In [94]:
locais_truncamento = [1, 10, 50, 100, 400]
n = len(locais_truncamento) + 1 # +1 para o original
ncols = 2
nrows = (n + 1) // 2 # número de linhas necessário
fig, axs = plt.subplots(nrows, ncols, figsize=(14, 3 * nrows))
axs = axs.flatten() # transforma em array 1D para indexação sequencial
# 1. Espectrograma original
plotar_espectrograma(espectrograma_db, taxa_amostragem, title="Original", ax=axs[0])
axs[0].set_title("Original (completa)")
# 2. Espectrogramas truncados por SVD
for i, k in enumerate(locais_truncamento):
idx = i + 1 # o índice 0 já foi usado para o original
U, S, Vt = decomposicaoSVD(espectrograma=espectrograma_db, k=k)
espectrograma_k = U @ S @ Vt
plotar_espectrograma(espectrograma_k, taxa_amostragem, title=f"Truncado (k = {k})", ax=axs[idx])
axs[idx].set_title(f"Truncado (k = {k})")
# Limpar eixos extras se houver
for j in range(len(locais_truncamento) + 1, len(axs)):
axs[j].axis('off')
plt.tight_layout()
plt.show()
4. Reconstrução do Espectrograma e do Áudio¶
In [231]:
# Número de valores singulares a serem mantidos na decomposição SVD.
k = np.where(valores_singulares == np.quantile(valores_singulares, 0.25))[0][0]
U, S, Vt = decomposicaoSVD(espectrograma_db, k)
espectrograma_resultante_db = U @ S @ Vt
# Convertendo de dB para magnitude linear
magnitude_resultante = 10 ** (espectrograma_resultante_db / 20)
# Reconstruindo a matriz STFT complexa
matriz_stft_resultante = magnitude_resultante * np.exp(1j * fase)
# Inverso da Transformada de Fourier de Curto Termo
_, audio_recuperado = istft(matriz_stft_resultante, fs=taxa_amostragem, nperseg=1024)
ipd.Audio(audio_recuperado, rate=sr)
Out[231]:
5. Compressão vs Qualidade — Análise Quantitativa¶
In [250]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# --- Parâmetros ---
valores_k = [5, 10, 25, 50, 100, 150, 200, 300, 400, 450, 500, 513]
caminho_audio = "audio/Confutatis.wav"
# --- Matriz original ---
espectrograma_db
linhas, colunas = espectrograma_db.shape
tamanho_original_elementos = linhas * colunas
tamanho_original_kb = espectrograma_db.nbytes / 1024
# --- Função de erro MSE ---
def erro_mse(original, reconstruido):
return np.mean((original - reconstruido) ** 2)
# --- Pré-cálculo da SVD completa ---
U, S, Vt = decomposicaoSVD(espectrograma_db)
# --- Inicialização de armazenamento ---
resultados = []
# --- Loop sobre os valores de k ---
for k in valores_k:
# Truncamento da SVD
U_k = U[:, :k]
S_k = S[:k, :k]
Vt_k = Vt[:k, :]
# Reconstrução do espectrograma_db
espectro_reconstruido = U_k @ S_k @ Vt_k
# Erro de reconstrução
erro = erro_mse(espectrograma_db, espectro_reconstruido)
# Tamanho em elementos (teórico)
tamanho_elementos = k * (linhas + colunas + 1)
taxa_elementos = tamanho_elementos / tamanho_original_elementos
# Tamanho real em KB (efetivo)
tamanho_kb = sum([matriz.nbytes for matriz in [U_k, S_k, Vt_k]]) / 1024
taxa_kb = tamanho_kb / tamanho_original_kb
# Armazenar resultado
resultados.append({
'k': k,
'Erro MSE': erro,
'Tamanho Comprimido (KB)': tamanho_kb,
'Taxa de Compressão (KB)': taxa_kb,
'Taxa de Compressão (Elementos)': taxa_elementos
})
# --- Criar DataFrame com resultados ---
df_resultados = pd.DataFrame(resultados)
# --- Exibir a tabela ordenada ---
display(df_resultados.sort_values(by='k'))
# --- Gráficos lado a lado ---
fig, axs = plt.subplots(1, 2, figsize=(15, 5))
# Gráfico 1: Erro MSE
axs[0].plot(df_resultados['k'], df_resultados['Erro MSE'], marker='o', color='#f4a261')
axs[0].set_title('Erro de Reconstrução (MSE) vs k', fontsize=12)
axs[0].set_xlabel('Número de Componentes k')
axs[0].set_ylabel('Erro Quadrático Médio (MSE)')
axs[0].grid(True, linestyle='--', alpha=0.5)
axs[0].set_xticks(valores_k)
axs[0].tick_params(axis='x', rotation=45)
# Gráfico 2: Taxa de Compressão (KB)
axs[1].plot(df_resultados['k'], df_resultados['Taxa de Compressão (KB)'], marker='o', color='#457b9d')
axs[1].axhline(1.0, linestyle='--', color='red', alpha=0.6, label='Tamanho Original')
axs[1].set_title('Taxa de Compressão Real (KB) vs k', fontsize=12)
axs[1].set_xlabel('Número de Componentes k')
axs[1].set_ylabel('Taxa de Compressão (comprimido / original)')
axs[1].grid(True, linestyle='--', alpha=0.5)
axs[1].legend()
axs[1].set_xticks(valores_k)
axs[1].tick_params(axis='x', rotation=45)
plt.suptitle('Análise de Compressão via SVD', fontsize=14, y=1.05)
plt.tight_layout()
plt.show()
| k | Erro MSE | Tamanho Comprimido (KB) | Taxa de Compressão (KB) | Taxa de Compressão (Elementos) | |
|---|---|---|---|---|---|
| 0 | 5 | 4.45e+01 | 146.07 | 0.01 | 0.01 |
| 1 | 10 | 3.96e+01 | 292.34 | 0.02 | 0.02 |
| 2 | 25 | 3.29e+01 | 732.32 | 0.05 | 0.05 |
| 3 | 50 | 2.80e+01 | 1469.53 | 0.11 | 0.10 |
| 4 | 100 | 2.15e+01 | 2958.59 | 0.21 | 0.21 |
| 5 | 150 | 1.65e+01 | 4467.19 | 0.32 | 0.31 |
| 6 | 200 | 1.24e+01 | 5995.31 | 0.43 | 0.42 |
| 7 | 300 | 6.47e+00 | 9110.16 | 0.65 | 0.63 |
| 8 | 400 | 2.71e+00 | 12303.12 | 0.88 | 0.84 |
| 9 | 450 | 1.35e+00 | 13928.91 | 1.00 | 0.94 |
| 10 | 500 | 2.34e-01 | 15574.22 | 1.12 | 1.05 |
| 11 | 513 | 1.24e-09 | 16005.20 | 1.15 | 1.07 |